Export as excel tables
for(ds in 1:2) {
library(openxlsx)
Up.wb <- createWorkbook()
for(i in 1:length(tbls)) {
res = filtered.tables[[i]][[ds]]
res = res[res$logFC > 0, ]
res = cbind(data.frame(Gene = rownames(res)), res)
res = res[order(res$t, decreasing = T), ]
n = names(filtered.tables)[[i]] #str_replace(arch.names[arch.order[i]], "/", "-")
addWorksheet(wb=Up.wb, sheetName = n)
writeData(Up.wb, sheet = n, res)
}
saveWorkbook(Up.wb, sprintf("~/PFC_v3/DE_genes_up_%s_subtype.xlsx", names(filtered.tables[[i]])[[ds]]), overwrite = TRUE)
library(openxlsx)
Down.wb <- createWorkbook()
for(i in 1:length(tbls)) {
res = filtered.tables[[i]][[ds]]
res = res[res$logFC < 0, ]
res = cbind(data.frame(Gene = rownames(res)), res)
res = res[order(res$t, decreasing = F), ]
n = names(filtered.tables)[[i]] #str_replace(arch.names[arch.order[i]], "/", "-")
addWorksheet(wb=Down.wb, sheetName = n)
writeData(Down.wb, sheet = n, res)
}
saveWorkbook(Down.wb, sprintf("~/PFC_v3/DE_genes_down_%s_subtype.xlsx", names(filtered.tables[[i]])[[ds]]), overwrite = TRUE)
}
combined.analysis.tables = lapply(names(filtered.tables), function(celltype) {
print(celltype)
tbls = filtered.tables[[celltype]]
gene.tbls = lapply(1:nrow(tbls[[1]]), function(i) {
dfs = lapply(1:length(tbls), function(k) tbls[[k]][i, ])
df = do.call("rbind", dfs)
})
names(gene.tbls) = tbls[[1]]$gene
combined.analysis.tbl = do.call(rbind, lapply(names(gene.tbls), function(gene){
x = suppressWarnings(metafor::rma(yi=logFC, sei=se, data = gene.tbls[[gene]], method="FE"))
combined.tbl = data.frame( gene = gene,
logFC = x$beta,
se = x$se,
tstat = x$zval,
P.Value = x$pval)
return(combined.tbl)
}))
rownames(combined.analysis.tbl) = names(gene.tbls)
combined.analysis.tbl = combined.analysis.tbl[order(combined.analysis.tbl$P.Value), ]
return(combined.analysis.tbl)
})
names(combined.analysis.tables) = names(filtered.tables)
DF = do.call(rbind, combined.analysis.tables)
DF$adj.P.Val = p.adjust(DF$P.Value, "fdr")
combined.analysis.tables = split(DF, unlist(lapply(names(combined.analysis.tables), function(celltype) rep(celltype, nrow(combined.analysis.tables[[celltype]])))))
readr::write_rds(combined.analysis.tables, "~/PFC_v3/meta_analysis_diff_results_subtypes.rds")
resDE = readr::read_rds("~/PFC_v3/filtered_resDE.rds")
filtered.tables = readr::read_rds("~/PFC_v3/individual_diff_results_filtered.rds")
combined.analysis.tables = readr::read_rds("~/PFC_v3/meta_analysis_diff_results.rds")
library(openxlsx)
Up.wb <- createWorkbook()
for(i in 1:length(combined.analysis.tables)) {
res = combined.analysis.tables[[i]]
res = res[(res$logFC > 0.1) & (res$P.Value <= 0.05), ]
res = res[order(res$t, decreasing = T), ]
n = names(combined.analysis.tables)[[i]] #str_replace(arch.names[arch.order[i]], "/", "-")
addWorksheet(wb=Up.wb, sheetName = n)
writeData(Up.wb, sheet = n, res)
}
saveWorkbook(Up.wb, sprintf("~/PFC_v3/DE_genes_up_%s_subtype.xlsx", "combined"), overwrite = TRUE)
library(openxlsx)
Down.wb <- createWorkbook()
for(i in 1:length(combined.analysis.tables)) {
res = combined.analysis.tables[[i]]
res = res[(res$logFC < -0.1) & (res$P.Value <= 0.05), ]
res = res[order(res$t, decreasing = F), ]
n = names(combined.analysis.tables)[[i]] #str_replace(arch.names[arch.order[i]], "/", "-")
addWorksheet(wb=Down.wb, sheetName = n)
writeData(Down.wb, sheet = n, res)
}
saveWorkbook(Down.wb, sprintf("~/PFC_v3/DE_genes_down_%s_subtype.xlsx", "combined"), overwrite = TRUE)
DE.sc = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
tstats = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
logFC = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
logPvals = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
rownames(DE.sc) = rownames(tstats) = rownames(logFC) = rownames(logPvals) = rownames(pb.logcounts)
colnames(DE.sc) = colnames(tstats) = colnames(logFC) = colnames(logPvals) = names(combined.analysis.tables)
limma_trend_mean.scores = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
Up.genes = vector("list", length(combined.analysis.tables))
Down.genes = vector("list", length(combined.analysis.tables))
rownames(limma_trend_mean.scores) = rownames(pb.logcounts)
names(Up.genes) = names(Down.genes) = colnames(limma_trend_mean.scores) = names(combined.analysis.tables)
for(i in 1:length(combined.analysis.tables)) {
print(i)
tbl = combined.analysis.tables[[i]]
tstats[tbl$gene, names(combined.analysis.tables)[[i]]] = tbl$tstat
logFC[tbl$gene, names(combined.analysis.tables)[[i]]] = tbl$logFC
logPvals[tbl$gene, names(combined.analysis.tables)[[i]]] = -log10(tbl$adj.P.Val)
DE.sc[tbl$gene, names(combined.analysis.tables)[[i]]] = tbl$tstat
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
limma_trend_mean.scores[is.na(limma_trend_mean.scores)] = 0
Up.genes = lapply(combined.analysis.tables, function(combined.analysis.tbl) {
combined.analysis.tbl$gene[(combined.analysis.tbl$logFC > 0.1) & (combined.analysis.tbl$adj.P.Val < 0.05)]
})
Down.genes = lapply(combined.analysis.tables, function(combined.analysis.tbl) {
combined.analysis.tbl$gene[(combined.analysis.tbl$logFC < -0.1) & (combined.analysis.tbl$adj.P.Val < 0.05)]
})
DE.new = list(DE.sc = DE.sc, tstats = tstats, logFC = logFC, logPvals = logPvals, Up.genes = Up.genes, Down.genes = Down.genes)
saveRDS(DE.new, "~/PFC_v3/DE_genes_pseudobulk_final_subtypes.rds")

ALS.DEs = read.table("~/als/MAGMA_DEG_input.tsv", sep = "\t")
suppressWarnings(ids <- AnnotationDbi::mapIds(org.Hs.eg.db, keys = ALS.DEs$V2, keytype = "ENSEMBL", column = "SYMBOL", multiVals = "first"))
'select()' returned 1:many mapping between keys and columns
ids[is.na(ids)] = ""
ids = as.character(ids)
ALS.genes = split(ids, ALS.DEs$V1)
X = do.call(cbind, lapply(ALS.genes, function(gs) as.numeric(rownames(variant_gene_scores) %in% gs)))
rownames(X) = rownames(variant_gene_scores)
# gg = apply(variant_gene_scores, 2, function(x) rownames(X)[scale(x) > 1])
EN = assess.geneset.enrichment.from.scores(variant_gene_scores, as(X, "sparseMatrix"))
logPvals = EN$logPvals
colnames(logPvals) = colnames(variant_gene_scores)
require(ComplexHeatmap)
pdf("~/als/MAGMA.pdf", height = 48)
Heatmap(logPvals)
dev.off()
null device
1
require(ComplexHeatmap)
pdf("~/als/MAGMA.pdf", height = 48)
ComplexHeatmap::Heatmap(Res[, -5], rect_gp = gpar(col = "black"))
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
closing unused connection 8 (~/als/magma/hmagmaAdultBrain__sz3_DE.gsa.out)
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
closing unused connection 7 (~/als/magma/hmagmaAdultBrain__ms_DE.gsa.out)
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
closing unused connection 6 (~/als/magma/hmagmaAdultBrain__alzKunkleNoapoe_DE.gsa.out)
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
closing unused connection 4 (~/als/magma/hmagmaAdultBrain__alz2noapoe_DE.gsa.out)
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
closing unused connection 3 (~/als/magma/hmagmaAdultBrain__als_DE.gsa.out)
dev.off()
null device
1

visualize.markers(In.ace, c("CALB1", "CALB2", "CUX2", "NOS1", "NPY"))
[1] "Markers Visualized: CALB1, CALB2, CUX2, NOS1, NPY"
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]





# Rosehip
# VIP
# REELIN
#
# PV-Basker
# PV-Chan
# SST

epi_markers$Sst_Calb1
[1] "HPGD+" "CALB1+" "FBN2+" "TAC3-" "TH-" "GXYLT2-" "ADGRG6-"
epi_annot = readr::read_rds("~/PFC_v3/epilepsy_paper_cell_annot.rds")
epi_markers = readr::read_rds("~/PFC_v3/epilepsy_paper_markers.rds")
In.cell.annot = annotate.cells.using.markers(In.ace, epi_markers[20:48])
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
cc = grep("Sst", colnames(In.cell.annot$Enrichment))
sapply(cc, function(i) {
plot(plot.ACTIONet.gradient(In.ace, x = In.cell.annot$Enrichment[, i], alpha_val = 0)+ggtitle(colnames(In.cell.annot$Enrichment)[[i]]))
})
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
data List,0 List,0 List,0 List,0 List,0 List,0 List,0 List,0 List,0 List,0
layers List,1 List,1 List,1 List,1 List,1 List,1 List,1 List,1 List,1 List,1
scales ? ? ? ? ? ? ? ? ? ?
mapping List,0 List,0 List,0 List,0 List,0 List,0 List,0 List,0 List,0 List,0
theme List,10 List,10 List,10 List,10 List,10 List,10 List,10 List,10 List,10 List,10
coordinates ? ? ? ? ? ? ? ? ? ?
facet ? ? ? ? ? ? ? ? ? ?
plot_env ? ? ? ? ? ? ? ? ? ?
labels List,6 List,6 List,6 List,6 List,6 List,6 List,6 List,6 List,6 List,6













annotate.profile.using.markers <- function (feature.scores, marker.genes, rand.sample.no = 1000)
{
if (is.matrix(marker.genes) | is.sparseMatrix(marker.genes)) {
marker.genes = apply(marker.genes, 2, function(x) rownames(marker.genes)[x >
0])
}
specificity.panel = feature.scores
GS.names = names(marker.genes)
if (is.null(GS.names)) {
GS.names = sapply(1:length(GS.names), function(i) sprintf("Celltype %s",
i))
}
markers.table = do.call(rbind, lapply(names(marker.genes),
function(celltype) {
genes = marker.genes[[celltype]]
if (length(genes) == 0) {
err = sprintf("No markers left.\n")
stop(err, call. = FALSE)
}
signed.count = sum(sapply(genes, function(gene) grepl("\\+$|-$",
gene)))
is.signed = signed.count > 0
if (!is.signed) {
df = data.frame(Gene = genes, Direction = +1,
Celltype = celltype, stringsAsFactors = FALSE)
}
else {
pos.genes = (as.character(sapply(genes[grepl("+",
genes, fixed = TRUE)], function(gene) stringr::str_replace(gene,
stringr::fixed("+"), ""))))
neg.genes = (as.character(sapply(genes[grepl("-",
genes, fixed = TRUE)], function(gene) stringr::str_replace(gene,
stringr::fixed("-"), ""))))
df = data.frame(Gene = c(pos.genes, neg.genes),
Direction = c(rep(+1, length(pos.genes)), rep(-1,
length(neg.genes))), Celltype = celltype,
stringsAsFactors = FALSE)
}
}))
markers.table = markers.table[markers.table$Gene %in% rownames(specificity.panel),
]
if (dim(markers.table)[1] == 0) {
err = sprintf("No markers left.\n")
stop(err, call. = FALSE)
}
specificity.panel = specificity.panel[markers.table$Gene, ]
IDX = split(1:dim(markers.table)[1], markers.table$Celltype)
print("Computing significance scores")
set.seed(0)
Z = sapply(IDX, function(idx) {
markers = (as.character(markers.table$Gene[idx]))
directions = markers.table$Direction[idx]
mask = markers %in% colnames(specificity.panel)
A = as.matrix(specificity.panel[markers[mask], ])
sgn = as.numeric(directions[mask])
stat = A %*% sgn
rand.stats = sapply(1:rand.sample.no, function(i) {
rand.samples = sample.int(dim(specificity.panel)[2],
sum(mask))
rand.A = as.matrix(specificity.panel[, rand.samples])
rand.stat = rand.A %*% sgn
})
cell.zscores = as.numeric((stat - apply(rand.stats, 1,
mean))/apply(rand.stats, 1, sd))
return(cell.zscores)
})
Z[is.na(Z)] = 0
Labels = colnames(Z)[apply(Z, 1, which.max)]
Labels.conf = apply(Z, 1, max)
names(Labels) = rownames(specificity.panel)
names(Labels.conf) = rownames(specificity.panel)
rownames(Z) = rownames(specificity.panel)
out = list(Label = Labels, Confidence = Labels.conf, Enrichment = Z)
return(out)
}

names(curatedMarkers_human$Brain$PFC$Velmeshev2019$marker.genes)
[1] "AST" "Endothelial" "IN-PV" "IN-SST" "IN-SV2C" "IN-VIP" "L2/3"
[8] "L4" "L5/6" "L5/6-CC" "Microglia" "Neu-NRGN" "Oligodendrocytes" "OPC"

Ast.DE = annotate.profile.using.markers(t(spec), mm)
[1] "Computing significance scores"
Heatmap(Ast.DE$Enrichment)
